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Feedback  Control  for  Aerodynamics 


R.  Chris  Camphouse*  Seddik  M.  Djouadi  t  and  James  H.  Myatt  + 


1.  Abstract 

The  two-dimensional  Burgers  equation  is  used  as  a  sur¬ 
rogate  for  the  governing  equations  to  test  order-reduction 
and  control  design  approaches.  This  scalar  equation  is  se¬ 
lected  because  it  has  a  nonlinearity  that  is  similar  to  the 
Navier-Stokes  equation,  but  it  can  be  accurately  simulated 
using  far  fewer  states.  However,  the  number  of  states  re¬ 
quired  is  still  well  above  that  for  which  a  controller  can  be 
designed  directly.  Two  approaches  for  order  reduction  are 
used.  In  both  approaches,  proper  orthogonal  decomposi¬ 
tion  (POD),  also  know  as  Karhunen-Loeve  decomposition 
or  principal  component  analysis,  is  used  with  Galerkin  pro¬ 
jection.  In  the  first  method,  the  traditional  POD  approach 
of  selecting  the  modes  to  be  retained  in  the  reduced-order 
model  is  based  on  the  energy  content  of  the  modes.  In  the 
second  method,  balanced  truncation  is  used  to  select  the 
appropriate  modes.  Both  approaches  capture  the  dynamics 
of  the  input-output  system  and  are  used  for  control  design. 

2.  Introduction 

Aerodynamic  flow  control  is  a  research  area  of  great  in¬ 
terest  to  the  United  States  Air  Force  and  the  fluid  mechan¬ 
ics  community  in  general.  Recent  advances  in  actuators, 
sensors,  simulation,  and  experimental  diagnostics  bring  ap¬ 
plications  such  as  suppression  of  acoustic  tones  in  cavi¬ 
ties,  separation  control  for  high  lift,  and  trajectory  control 
without  moving  hinged  surfaces  within  reach.  However, 
many  applications  require  the  integration  of  feedback  con¬ 
trol  because  of  the  need  for  robustness  to  flight  condition 
and  vehicle  attitude,  precision  tracking,  overcoming  low- 
fidelity  models,  or  moving  a  system  away  from  a  stable 
solution  or  limit  cycle  as  efficiently  as  possible.  Feedback 
control  strategies  in  which  the  bandwidth  of  the  controller 
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is  commensurate  with  the  time  scales  of  the  aerodynamics 
are  attractive  because  they  offer  the  possibility  of  improved 
performance  and  reduced  control  power  required  through 
control  of  unstable  structures  in  the  flow  field.  Unfortu¬ 
nately,  models  that  capture  the  relevant  dynamics  of  the 
input-output  system  and  are  amenable  to  control  design  are 
difficult  to  develop. 

The  governing  equations  for  a  compressible  fluid  are 
partial  differential  equations  -  1)  the  Navier-Stokes  equa¬ 
tion  for  momentum,  2)  the  continuity  equation,  and  3)  the 
energy  equation.  As  an  illustration  of  the  complexity  of 
the  governing  equations,  the  dimensionless  Navier-Stokes 
equation  for  momentum  [1]  for  an  incompressible,  New¬ 
tonian  fluid  with  a  few  simplifying  assumptions  such  as  no 
body  force  is 

^  +  (m- V)m  +  V/?=  -J-Aiq  (1) 

dt  Re 

where  u  =  u(t,x,y,z)  is  the  velocity,  p  —  p(t,x,y,z)  is  the 
pressure,  and  the  dimensionless  Reynolds  number  Re  is  a 
measure  of  the  ratio  of  inertia  forces  to  viscous  forces.  As 
Re  increases  with  vehicle  size  and  speed,  the  effect  of  the 
linear  terms  is  diminished  and  the  nonlinear  terms  become 
more  dominant.  This  greatly  increases  challenges  of  both 
the  modeling  and  control  problems  for  vehicles  at  realistic 
flight  conditions. 

Computational  fluid  dynamics  (CFD)  simulations  can 
provide  solutions  to  a  discretized  form  of  the  Navier- 
Stokes.  However,  accurate  simulations  for  simple  shapes 
such  as  two-dimensional  airfoils  can  require  several  thou¬ 
sand  states  and  therefore  are  not  directly  useful  for  con¬ 
trol  design  due  to  the  extremely  high  order  of  the  system. 
Simulations  for  a  full  vehicle  can  require  over  one  million 
states.  The  large  number  of  states  is  necessary  to  capture 
important  flow  features  that  occur  at  extremely  small  spa¬ 
tial  scales.  Another  challenge  for  designing  control  laws 
for  flow  control  is  that  the  Navier-Stokes  equations  in  tra¬ 
ditional  form  are  not  affine.  Therefore,  it  is  necessary  to 
separate  the  portions  of  the  discretized  system  where  the 
control  input  enters  the  system  for  control  design.  In  ad¬ 
dition,  the  location  of  the  actuator  is  often  limited  to  the 
boundary.  This  eliminates  the  possibility  of  full-state  feed¬ 
back. 
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3.  Model  Problem 

3.1.  Motivation:  Boundary  Feedback  Control 

Reduced  order  modelling  has  received  significant  atten¬ 
tion  by  the  research  community  in  recent  years.  Proper  or¬ 
thogonal  decomposition  (POD)  has  been  investigated  ex¬ 
tensively  [2]- [10]  as  a  potential  technique  due  to  its  signif¬ 
icant  order  reduction  capability.  Full  order  system  model 
behavior  described  by  thousands  of  states  can  often  be  cap¬ 
tured  with  a  POD  model  composed  of  dozens  of  states  or 
less.  The  implications  for  control  law  design  are  obvious. 
Applications  requiring  exceedingly  large  systems  for  accu¬ 
rate  simulation  result  in  intractable  feedback  control  prob¬ 
lems.  A  particular  example  is  feedback  control  of  aerody¬ 
namic  fluid  flows.  Accurate  simulations  done  with  com¬ 
putational  fluid  dynamics  typically  require  discretized  flow 
models  describing  thousands  of  states,  usually  millions  in 
the  case  of  three-dimensional  turbulent  flow.  The  system¬ 
atic  development  of  feedback  control  laws  from  systems  of 
such  large  dimension  is  currently  an  intractable  problem. 
Reduction  of  system  order  must  be  done  if  feedback  con¬ 
trol  law  design  for  these  systems  is  to  be  feasible. 

In  many  practical  applications,  boundary  actuation  is  a 
requirement.  For  example,  control  of  flow  separation  over 
an  airfoil  requires  that  actuation  and  sensing  be  done  on  the 
airfoil  surface.  The  possibility  of  unmodelled  dynamics  in 
the  system,  or  dynamics  lost  in  the  order  reduction  process, 
make  feedback  control  a  requirement.  Systematic  develop¬ 
ment  of  boundary  feedback  control  laws  from  POD  models 
has  remained  an  elusive  problem.  Controls  are  often  in¬ 
cluded  in  the  reduced  model  in  an  ad-hoc  way,  specified  to 
be  distributed  over  a  subregion  of  the  domain  interior,  or 
simply  specified  through  open-loop  forcing. 

In  this  paper,  we  utilize  the  weak  formulation  of  reduced 
order  models  obtained  via  POD  and  Galerkin  projection. 
Use  of  the  weak  form  permits  separation  of  the  boundary 
input  in  the  reduced  model,  allowing  the  boundary  control 
to  enter  the  reduced  model  equations  explicitly.  We  present 
these  concepts  using  a  nonlinear  convective  system  defined 
over  an  obstacle  geometry.  The  resulting  model  problem 
captures  many  of  the  difficulties  associated  with  feedback 
control  of  fluid  flows  with  a  greatly  reduced  computational 
workload. 

3.2.  Distributed  Parameter  System 

Let  C  R2  be  the  rectangle  given  by  (a,  b]  x  (c,  d).  Let 
0,2  C  be  the  rectangle  given  by  [01,02]  x  [^1,^2]  where 
a  <  a  1  <  02  <  b  and  c  <  b\  <b2<d.  The  problem  domain, 
O,  is  given  by  O  =  Qi  \  O2.  In  this  configuration,  O2  is  the 
obstacle.  Dirichlet  boundary  controls  are  located  on  the  ob¬ 
stacle  bottom  and  top,  denoted  by  Tg  and  Ft,  respectively. 


Figure  1 .  Problem  Geometry. 


The  dynamics  of  the  system  are  described  by  the  two- 
dimensional  Burgers  equation 


^w(t,x,y)  +  W  -F(w)  =  La  w{t,x,y) 
for  t  >  0  and  (x,y)  €  £!.  In  (2),  F(w)  has  the  form 


F(yv)  = 


Ci 


w2(t,x,y) 


C2 


w2(t,x,y) 


(2) 

(3) 


where  C\,  C2  are  nonnegative  constants.  This  system  has 
a  convective  nonlinearity  like  that  found  in  the  Navier- 
Stokes  partial  differential  equations  modeling  fluid  flows. 
The  quantity  Re ,  a  nonnegative  constant,  is  analogous  to 
the  Reynolds  number  in  the  Navier- Stokes  equations. 

To  complete  the  model  of  the  system,  boundary  condi¬ 
tions  must  be  specified  as  well  as  an  initial  condition.  For 
simplicity,  boundary  controls  are  assumed  to  be  separable. 
With  this  assumption,  we  specify  conditions  on  the  obsta¬ 
cle  bottom  and  top  of  the  form 


w(t,  Fb)  =  (4) 

w(t,FT)  =  uT(t)'FT{x).  (5) 

In  (4)-(5),  usit)  and  uj(t)  are  the  controls  on  the  bottom 
and  top  of  the  obstacle,  respectively.  The  profile  functions 
(jc)  and  TV(x)  describe  the  spatial  influence  of  the  con¬ 
trols  on  the  boundary.  A  parabolic  inflow  condition  is  spec¬ 
ified  of  the  form 


w(t,Fin)=f(y).  (6) 

At  the  outflow,  a  Neumann  condition  is  specified  according 

to 

-^w(t,rout)  =  o.  (7) 

For  notational  convenience,  denote  the  remaining  boundary 
as  Tf/.  We  require  that  values  be  fixed  at  zero  along  Tf/  as 
time  evolves.  The  resulting  boundary  condition  is  of  the 
form 

w(t,rv)  =  o.  (8) 

The  initial  condition  of  the  system  is  given  by 

w(0,x,y)  =  w0{x,y)  6  L2(Q.).  (9) 
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3.3.  POD  Basis  Construction 

The  snapshot  method  [1 1]  is  used  to  construct  a  low  or¬ 
der  POD  basis  for  the  distributed  parameter  system.  An 
ensemble  of  solution  snapshots  {5/(x)}^1  for  system  (2), 
(4)-(9)  is  generated  by  numerical  simulation.  In  (3),  we  set 
C\  =  1  and  C2  =  0  in  order  to  obtain  solutions  that  con- 
vect  from  left  to  right  for  positive  inflow  condition  f(y).  In 
addition,  we  specify  that  Re  =  300.  Finite-difference  spa¬ 
tial  discretization  is  done  as  in  [12].  A  uniform  Cartesian 
grid  is  constructed  with  spatial  step-size  h.  The  resulting 
discretized  system  describes  roughly  2,000  states.  As  it  is 
desired  that  the  POD  basis  spans  dynamics  introduced  by 
a  time- varying  boundary  control  input,  nontrivial  boundary 
conditions  are  specified  during  ensemble  generation.  In¬ 
puts  specified  are  of  the  form 

UB(t)  =  D s in (0. 25 12  25 )  uj(t)  =  0,  (10) 

us(t)  =  0  uj(t)  =  Dsin(0.25t2'25),  (11) 

for  D  =  —3,  -2,  -1.  The  steady  solution  arising  from  a  pos¬ 
itive  parabolic  inflow  condition  is  taken  as  the  initial  con¬ 
dition  in  each  simulation.  For  each  case  of  control  input 
listed  in  (10)  -  (11),  snapshots  are  taken  in  increments  of 
At  =  0.1  starting  from  t  —  0  and  ending  at  T  =  15.  The 
snapshots  resulting  from  each  case  are  combined  into  an 
overall  snapshot  set.  The  resulting  ensemble  consists  of 
roughly  900  snapshots. 

With  the  snapshot  ensemble  in  hand,  the  N  x  N  correla¬ 
tion  matrix  L  defined  by 

Lij  =  (Si,Sj )  (12) 

is  constructed.  In  this  work,  we  utilize  the  standard  L 2  (£2) 
inner  product 

(Si,sj)  =  Jsifjdz,  (13) 

where  S*  denotes  the  complex  conjugate  of  Sj,  in  the  con¬ 
struction  of  L. 

With  M  denoting  the  number  of  POD  modes  to  be 
constructed,  the  first  M  eigenvalues  of  largest  magnitude, 
of  L  are  found.  They  are  sorted  in  descending  or¬ 
der,  and  their  corresponding  eigenvectors  {v/}^1  are  cal¬ 
culated.  Each  eigenvector  is  normalized  so  that 

)N:-T  (14) 

The  orthonormal  POD  basis  set  {0/(x)}^1  is  con¬ 

structed  according  to 


N 


<MX)  =  L  Vij57(x), 

.7=1 


(15) 


where  vy  is  the  jth  component  of  v/. 

With  a  POD  basis  in  hand,  the  solution  w(t,x)  of  the  dis¬ 
tributed  parameter  model  is  approximated  as  a  linear  com¬ 
bination  of  POD  modes,  i.e., 

M 

w(t,x) «  £a,V)0,-(x).  (16) 

i=  1 

3.4.  Weak  Galerkin  Model 


We  now  develop  a  reduced  order  model  for  the  system 
described  by  (2),  (4)-(9).  Using  the  weak  formulation  of 
the  governing  equation  allows  us  to  extract  boundary  con¬ 
dition  information  prior  to  Galerkin  projection.  Galerkin 
projection  of  the  weak  system  onto  the  POD  basis  results 
in  a  system  of  ordinary  differential  equations  for  the  tem¬ 
poral  coefficients  with  explicit  control  input. 

Taking  the  inner  product  of  both  sides  of  (2)  with  the  i-th 
POD  mode  0/(x,y)  and  utilizing  Green’s  identities  results 
in  the  weak  formulation 


C  d 

Jn  diw^t,x,y^x,y^dK  = 

[  F  -V(j>i(x,y)dx-  f  (F(w)-n)<l>i(x,y)dA(x) 

JQ  JdCl 


1 

+  Re 


/  dd 


(Vw(t,x,y)  -n)<lH(x,y)dA(x) 


(17) 


1 

Re 


J  Vw(t,x,y)  ■V0i(x,y) ds 


where  n  denotes  the  unit  outward  normal. 

As  seen  in  (15),  each  POD  mode  is  a  linear  combination 
of  solution  snapshots.  From  (8),  snapshot  values  along  Tf/ 
are  specified  to  be  zero.  As  a  result,  POD  modes  are  zero 
along  Tf/.  Thus,  the  second  boundary  integral  in  (17)  is 
decomposed  as 


/  (Vw(t,x,y)  ■  n )<j>i(x,y)dA(x)  = 

Jan 

r“2  ( d  d  \ 

J  l  jdx 

fd  d 

-J  -^w(t,a,y)(j)i(a,y)dy, 


(18) 


where  condition  (7)  has  been  used  to  specify  that 

Ic  ^W^’b,y^^b’y^dy  =  0-  (19) 

In  a  similar  fashion,  the  remaining  boundary  integral  in 
(17)  is  decomposed  as 

[  (F  (w)  •  n)  0,-  (x,  y)dA  (x)  = 

JdQ. 

^  (w{t,b,y)2<l>i(b,y)  -  f(y)2<j>i(a, 


(20) 

y))  dy, 
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where  (6)  has  been  used  to  incorporate  the  inflow  condition 

m 

Control  inputs  and  the  Dirichlet  inflow  condition  are  not 
explicit  in  (18).  They  can  be  made  explicit  by  approximat¬ 
ing  partial  derivatives  along  the  boundary.  For  h  >  0,  we 
see  that 


d  ,  , 

2~w{t,x,bi)  K 

^  uB(t)x¥B(x)-w(t,x,bi-h) 

(21) 

h 

d 

-jj-w(t,x,b2)  s 

w(t,x,b2  +  h)  —  UT(t)x¥r(x) 

(22) 

h 

Txw{t’a'y)* 

„  w{t,a  +  h,y)-f(y) 

(23) 

h 

These  expressions  are  substituted  into  (18).  Approxi¬ 
mating  w(t,x,y)  as  a  linear  combination  of  POD  modes 
in  (17),  (18),  and  (20)  results  in  a  reduced  order  system 
model.  By  defining  fi  as 


1 

hRe 


(24) 


the  reduced  order  model  is  of  the  form 


ex — A.oc  Bu  A(oj)  F:  (25) 


N(a)  = 


1 

2 


fa  (EjL i  a;4>j(x’y))  £Mx,y) dx 
fc  (ijl  1  <Xj<l>j(b,y))  01  (b,y)dy 


1 

2 


//  (e"  1  aj<l>j(b,y))  <t>M(b,y)dy 


(28) 


Scinffy)  +  \f{y)2)<i>\{a,y)dy 


F  = 


(29) 


UcW(y)  +  {fiyf^M^a^dy \ 


Mx  1 


Projecting  the  initial  condition  wo(x,y)  onto  the  POD 
basis  results  in  an  initial  condition  for  the  reduced  order 
model  of  the  form 


a(  0)  =  do¬ 


do) 


3.5.  Model  Validation 


where 


MiJ)  = 

ra  2 


ra  2 

At[  /  (<l>j(x,b l  - h)<l>i(x,b1)  +  <j>j(x,b2  +  h)<j>i(x,b2))dx 

Ja\ 

-  f  <t>j(a  +  h,y)<pi(a,y)dy  +  h  J^V<j)j(x,y)  -V<l)j(x,y)]dx, 

(26) 


B  = 


Q  Mx,bl)^B(x)dx  Q  Mx,b2)^r(x)dx 


A* 


f^Mx,bi)^B(x)dx 


-I  Mx 2 

(27) 


Before  using  the  reduced  order  model  in  (25)  to  design 
feedback  control  laws,  we  first  verify  agreement  between 
the  reduced  and  full  order  systems.  We  utilize  the  ratio 


100 


(31) 


to  determine  how  many  POD  modes  to  include  in  the  re¬ 
duced  model.  The  POD  basis  is  optimal  in  an  energy  sense 
[13].  It  captures  the  mean  square  energy  of  the  snapshot 
ensemble  better  than  any  other  basis.  The  quantity  in  (31) 
provides  a  measure  of  the  ensemble  energy  that  is  cap¬ 
tured  by  the  POD  basis.  By  requiring  that  99.9%  of  the 
energy  contained  in  the  snapshot  ensemble  be  contained  in 
the  POD  basis,  we  calculate  the  smallest  value  of  M  such 
that  the  quantity  in  (31)  is  greater  than  or  equal  to  99.9.  For 
our  snapshot  ensemble,  the  smallest  value  of  M  satisfying 
this  relationship  is  15.  Thus,  we  include  15  POD  modes 
in  the  construction  of  the  reduced  model  in  (25).  The  first 
nine  modes  are  shown  in  Figure  2. 

We  now  compare  the  solution  obtained  from  the  re¬ 
duced  and  full  order  models  using  different  boundary  in¬ 
puts  that  those  used  during  ensemble  creation.  Boundary 
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Figure  2.  9  POD  Modes. 

inputs  specified  are  of  the  form 

uB(t)  =min(^~,  l)  ,  (32) 


We  first  verify  boundary  condition  agreement  between  the 
full  order  system  and  the  linear  combination  of  POD  modes 
given  by  (16)  with  the  modes  restricted  to  the  boundary. 
By  specifying  characteristic  functions  for  the  control  pro¬ 
file  functions  ¥#(*)  and  'Pj’(jc)  in  (4)-(5),  we  see  that 


M 


X>(f)fc(riO* 

i=  1 

aw(t,rB)  =  uB(t), 

(34) 

M 

I>(0fc(rY)s 

i=  1 

-  w(t,rT)  =  uT(t), 

(35) 

We  construct  the  linear  combinations  given  on  the  left  of 
(34)-(35)  and  compare  them  to  the  boundary  inputs  given 
by  (32)-(33).  The  results  are  shown  in  Figure  3.  In  Fig¬ 
ure  3,  dashed  lines  denote  the  linear  combination  of  POD 
modes  restricted  to  the  boundary.  Solid  lines  denote  the 
boundary  inputs  defined  by  (32)-(33).  As  can  be  seen  in 
Figure  3,  there  is  very  good  agreement  between  the  bound¬ 
ary  conditions  specified  for  the  full  order  system  and  the 
linear  combination  of  POD  modes  restricted  to  the  bound¬ 
ary.  To  further  validate  the  reduced  model,  we  project  the 
solution  from  the  full  order  simulation  at  each  time  step 
onto  the  POD  basis.  The  resulting  temporal  coefficients  are 
compared  to  those  predicted  from  the  reduced  order  model. 
The  results  obtained  for  the  first  five  temporal  coefficients 
are  shown  in  Figure  4.  In  that  figure,  solid  curves  denote 
values  of  temporal  coefficients  obtained  from  the  projec¬ 
tion.  Dashed  curves  denote  the  solution  of  the  reduced  or¬ 
der  model.  As  seen  in  Figure  4,  very  good  agreement  is 


Figure  3.  Boundary  Condition  Accuracy. 


Figure  4.  Projected  And  POD  Model  Coefficients. 

seen  between  the  full  and  reduced  solutions  even  though 
the  open-loop  input  considered  was  not  specifically  incor¬ 
porated  in  the  snapshot  ensemble. 

3.6.  Linear  Quadratic  Control 

The  system  given  by  (25),  (30)  is  linearized  about  the 
origin  yielding  a  state- space  equation  of  the  form 

a(t)  =Aa  +  Bu,  (36) 

a(0)  =  Ob-  (37) 

We  consider  the  tracking  problem  for  (36)-(37).  A  fixed 
reference  signal  wref(x)  is  specified  for  the  full  order  sys¬ 
tem.  Projecting  wr£?/(x)  onto  the  POD  basis  yields  tracking 
coefficients  for  the  reduced  order  model,  denoted  by  aref. 


5 


Reference  Function 


The  dynamics  of  the  linearized  model  under  tracking  con¬ 
trol  are  given  by 


a 

A 

0 

a 

B 

CCref 

— 

0 

0 

CCref 

+ 

0 

—  AX  T  Bu , 

where  we  have  defined  the  augmented  state  X  as 


X(t)  = 


a{t) 

OCref 


with  Xq  = 


cco 

CCref 


(38) 

(39) 

(40) 


To  formulate  the  control  problem,  we  consider  the  y- 
shifted  linear  quadratic  regulator  (LQR)  cost  function 


J(ao,u)  =  J  {( a  -  aref)T  Q(a  -  aref)  +  uT  Ru}  e2yt  dt . 

(41) 

In  (41),  Q  is  a  diagonal,  symmetric,  positive  semi-definite 
matrix  of  state  weights.  R  is  a  diagonal,  symmetric,  pos¬ 
itive  definite  matrix  of  control  weights.  The  quantity  y,  a 
nonnegative  constant,  is  an  additional  parameter  that  pro¬ 
vides  added  robustness  in  the  control  [15].  The  optimal 
control  problem  we  consider  is  to  minimize  (41)  over  all 
controls  u  €  L2( 0,°o)  subject  to  the  constraints  (38)-(40). 

For  a  controllable  system,  the  LQR  problem  has  a 
unique  solution  of  the  form 


uopt  =  -KX  (42) 

=  -\Kx  K2\X  (43) 

=  -[/?_1firnn  r~1bt nl2]x,  (44) 

where  Tin  is  the  unique  symmetric,  non-negative  solution 

of  the  algebraic  Riccati  equation 

(a + y/)rnu  +  nu  (a + y/)  -  n„B/?-1srn1i  +  q  =  o. 

(45) 

The  matrix  H\2  in  (44)  satisfies  the  equation 

[(A  +  ylf  -nnBR-1BT]nu  =  Q-  (46) 

The  feedback  control  obtained  from  the  linearized 
model  is  placed  into  the  nonlinear  state-space  equation.  As 
discussed  in  [14]- [15],  the  resulting  closed-loop  nonlinear 
system  is  of  the  form 

X  =  (A  —  BK)X  +  [N(a)  0]r  +  [F  0]r,  (47) 

*(0)=X0.  (48) 


3.7.  Closed-Loop  Results 

The  tracking  LQR  problem  requires  the  specification 
of  the  reference  signal  wref(x).  In  the  results  that  follow, 
wref(x)  is  defined  as  the  unactuated  steady  solution  for  the 
case  Re  =  50.  To  obtain  the  reference  function  for  the 


Figure  5.  Reference  Signal. 


reduced  model,  we  project  wref(x)  onto  the  fifteen  POD 
modes.  The  values  obtained  are  used  as  tracking  coeffi¬ 
cients  in  the  reduced  order  control  problem.  The  reference 
signal  obtained  by  projecting  wref(x)  onto  the  POD  basis 
is  shown  in  Figure  5. 

The  values  specified  in  the  control  formulation  are  Q  = 
3500 /15x15,  R  =  I2x2,  and  y=  0.41.  We  specify  wo(x)  to  be 
the  steady-state  solution  for  the  case  of  Re  =  300.  There¬ 
fore,  the  solution  will  remain  at  wo(x)  until  the  boundary 
condition  at  the  actuator  location  is  altered  by  the  con¬ 
troller.  The  controlled  reduced  order  model  is  shown  in 
Figure  6.  By  comparing  the  controlled  solution  of  Figure  6 


t  *  0.05  I  =  2.S5 


Figure  6.  Controlled  POD  Model. 


to  the  reference  function  in  Figure  5,  it  is  apparent  that  the 
reduced  order  control  satisfies  the  control  objective  quite 
well.  Significant  tracking  is  achieved. 
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4.  Balanced  Truncation  and  H°°  Control 


where 


4.1.  Balanced  Truncation 

Balanced  truncation  is  a  simple  and  popular  model  re¬ 
duction  technique,  which  can  be  described  as  follows  [20, 
21,  22,  23]:  Suppose  we  have  a  stable  linear  time  invari¬ 
ant  (LTI)  system  described  by  the  following  n-dimensional 
state  space  equation 

x(t)  =  Ax(t)  +  Bu(t)  (49) 

y(t)  =  Cx(t) 

where  x(t)  is  the  n  x  1 -state  vector  of  the  system,  u(t)  is  an 
m  x  1 -input  vector,  and  y(t)  is  an  p  x  1 -output  or  measure¬ 
ment  vector.  A,  B ,  and  C  are  constant  matrices  of  appropri¬ 
ate  dimensions. 

The  underlying  idea  of  balanced  truncation  is  to  take  into 
account  both  the  input  and  output  signals  of  the  system 
when  deciding  which  states  to  truncate  with  appropriate 
scaling.  The  latter  is  performed  by  transforming  the  con¬ 
trollability  and  observability  gramians,  denoted  Wc  and  W0 
respectively,  so  that  they  are  equal  and  diagonal. 

The  controllability  and  observability  gramians  satisfy  the 
following  Lypaunov  equations  [21] 

AWc  +  WcAt  +BBt  =  0  (50) 

AtW0  +  W0A+CtC  =  0  (51) 

The  controllability  and  observability  gramians  can  be  rep¬ 
resented  as 

Wc  =  J  eA‘BBTeATtdt 

W0  =  [  eAl,CTCeA,dt  (52) 

Jo 

Computing  a  state  balancing  transformation  M  is  achieved 
by  first  calculating  the  matrix  [21,  22] 

Wco  =  WcWo  (53) 

and  determining  its  eigenmodes 

Wco=MAM~1  (54) 

Let 

z(t):=M~lx(t)  (55) 

then  the  resulting  transformed  state  space  is 

z(t)  =  Az(t)+Bu(t)  (56) 

y(t)  =  Cz(t ) 


A  :=  M~lAM 
B  :=  M~lB 
C  :=  CM 

The  transformation  M  is  chosen  such  that  the  controllabil¬ 
ity  and  observability  gramians  for  the  transformed  system 
satisfy 

WC=W0=M~1WCM~1T  =MtW0M=:I,  (57) 

where  E  is  a  diagonal  matrix  that  satisfies  E2  =  A,  and  the 
diagonal  elements  of  £,  oj  s,  are  known  as  the  Hankel  sin¬ 
gular  vales  of  the  system,  i.e., 

E  =  diag{<7i,  <72,  ,  on}  (58) 

where  (7/  are  the  Hankel  singular  values  of  the  system  G 
arranged  in  non-increasing  order 

G\  >  02  >  •  >  cn  >  0  (59) 

In  balanced  truncation  only  states  corresponding  to  large 
Hankel  singular  values  are  retained.  Small  Hankel  singular 
values  correspond  to  states  which  are  deemed  weakly  con¬ 
trollable  and  weakly  observable,  and  therefore  deleted  from 
the  state-space  model.  For  instance,  if  the  first  nr  states  are 
retained  then  the  resulting  transformation  is  given  by 

Mr  =  [Ir  0 \M  (60) 

where  Ir  is  the  nr  x  nr  identity  matrix.  The  reduced  order 
model  is  obtained  by  letting 

Xr=\lr  0 ]MX  (61) 


as  follows 

xr(t)  = 

Vr 

0  ]M~lAM  o 

Xr{t)  +  [Ir  0  }M~lBu(t) 

yr(t)  = 

CM 

0  xr(0 

(62) 

Let 

Ar 

1 

o 

II 

lAM  Q 

By 

II 

i 

lB 

Cr 

:=  CM  J 

(63) 

The  error  bound  for  the  output  is  given  by 

lbW->v(0ll2<2  Y,  °ilK?)l|2>  Vm  €  L2  (64) 

nr+ 1 
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Balanced  truncation  is  optimal  in  a  precise  sense  [26].  To 
see  this  define  a  causal  bounded  input-output  operator  G 
acting  on  L2  (— oo5  oo)  into  L2  (— oo,  oo)  described  by  the  con¬ 
volution  [22,  23] 

(< Gu)(t):=  J'  CeA(-,~^Bu(x)dx  (65) 

Now,  define  the  Hankel  operator 

rG:  L2(- oo,  0]^L2[0,  oo) 


More  explicitly  [26], 

r Gnru(t)=l  Crekri-t-x)Bru{x)dx,  for  f  >  0  (73) 

and  the  optimal  index  is  given  by 

av  =  l|rG-rGj|Hs  (74) 

Optimality  of  balanced  truncation  seems  to  be  missing 
in  the  literature.  In  fact,  it  has  been  widely  claimed  that 
balanced  truncation  is  not  optimal  in  any  sense  [20,  23,  24]. 


of  G  by 

rG:=P+G\L2{_^0]  (66) 

where  G\l2^_oo  0j  denotes  the  restriction  of  G  to  L2(— oo,  0], 

and  P+  is  the  orthogonal  projection  acting  from  L2  (— oo,  oo) 
into  L2[ 0,  oo),  i.e.,  P+  is  the  truncation  operator 

P+m  =  {  f0{t)  ,  fit)  £  L2(-oo ,  oo)  (67) 

Then,  the  Hankel  operator  Tg  can  be  written  as 

r0 

Tcu(t )  =  I  Ce^^~x^ Bu(z)dz ,  for  t  >  0  (68) 

The  Hankel  operator  Tg  maps  past  inputs  to  future  outputs. 
Expression  (68)  shows  that  the  Hankel  operator  Tg  is  an 
integral  operator  mapping  L2(— oo,  0)  into  L2[ 0,  oo),  with 
kernel  the  impulse  response  k(t,  t)  defined  by 

k(t,  x)  :=  CeM,~T>B,  x  <  0,  t  >  0  (69) 

The  Hankel  operator  Tg  has  finite  rank  k  <  n,  that  is,  its 
range  has  finite  dimension  k  <  n  [21,  23],  and  therefore  be¬ 
longs  to  the  Hilbert- Schmidt  class  of  operators  acting  from 
L2(— oo,  0]  intoL2[0,  oo)  [25].  Its  Hilbert- Schmidt  norm  is 
defined  as 


In  terms  of  kernel  approximation,  balanced  trunca¬ 
tion  is  a  particular  case  of  POD  in  the  sense  that  the  kernel 
we  want  to  approximate  is  the  impulse  response  of  the 
system  k(t ,  t)  defined  in  (71).  The  optimization  index  iinr 
can  then  be  written  as  in  POD  [26] 

Uoo  m  |  nr  1 2 

J  \k(t,  T)-  L-/iMSi(T)|  dxdt 

:  fi  €  L2[0,  oo);  g[  £  L2{—°°,  0])  (75) 

=  J  J°  \k(t,  x)-CreAr{t~'c)Br^dxdt  (76) 

Expressions  (74)  and  (76)  show  that  balanced  truncation 
is  optimal  in  the  sense  of  optimal  approximation  in  the 
Hilbert- Schmidt  norm  of  the  Hankel  operator  Tg,  and 
optimal  in  the  sense  of  the  ||  •  ||  2 -norm  of  kernels  corre¬ 
sponding  to  impulse  responses  of  linear  time-invariant 
systems  defined  over  [0,  00)  x  (—00,  0]. 

The  linear  time-invariant  system  framework  allows 
the  exact  computations  of  the  optimal  lower  order  model 
approximation.  This  contrasts  with  POD  which  uses 
simulation  data  and  particular  open-loop  inputs  to  generate 
snapshots. 

4.2.  Application  to  the  Weak  Galerkin  Model 


llrdlns  =  jTjf  \k{t,x)\2dxdt  (70) 

Next,  consider  the  optimal  distance  minimization  de¬ 

fined  in  (71),  which  consists  of  optimally  approximating  in 
the  Hilbert- Schmidt  norm  the  Hankel  operator  Tg  by  an¬ 
other  Hankel  operator  Tcnr  of  lesser  rank,  say  nr  <  k ,  in 
other  words 

Pnr  :=nfin||rG-rG  ||hs  (71) 

nr<k 

It  turns  out  that  the  minimizer  in  (71)  is  the  Hankel  oper¬ 
ator  with  kernel  the  impulse  response  of  the  reduced  order 
model  (62)  [26],  i.e., 

Cr^'-^Br  (72) 


Our  approach  is  to  construct  an  approximately  balanced 
realization  to  (25).  This  is  carried  out  by  first  linearizing 
(25)  around  a®.  The  state  space  and  output  equations  have 
the  form 

a(t)  =  Aoc(t)+Bu(t ),  a(0)  =  Oo  (77) 
y(t)  =  a(t)  (78) 


where  A  is  a  matrix  having  the  same  dimension  as  A,  and  is 
given  by 


A 

B 


d(Aa  +  Bu(t)  +N(a)+F ) 
da 

a=ao 

d  (Aa  -\-Bu(t)  -\-N(a)  +F)  ^ 

~  Ihi  =t 


(79) 

(80) 
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In  this  model  the  dimension  of  the  state  vector  a  is  40 
which  corresponds  to  40  POD  modes.  A  balanced  realiza¬ 
tion  is  first  computed  and  27  states  truncated,  i.e.,  only  the 
states  corresponding  to  the  13  largest  Hankel  singular  are 
kept  in  the  model.  This  results  in  a  13 -dimensional  state- 
space  model 

i(t)  —  Az(t)+Bu(t)  (81) 

y(t)  =  cz(t ) 

where  z(t)  €  R.13,  B  6  R13x2,  and  C  6  R40x13. 

The  first  8  POD  modes  corresponding  are  shown  in 
Figure  7.  We  project  the  solution  from  the  full  order 


The  top  control 


The  bottom  control 
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0  0.Z  0.4  0.6  o.a 

Mode  3 


i 


0  0.2  0.4  0.6  o.a 

Mode  5 


Mode  2 
0.4 

0.2  : 

0 - 

0  0.2  0  4  0.6  0  6 

Mode  4 

0.4  ^0^^ 

0.2  _ 

0 . 

o  o.2  0  4  o.e  oa 

Mode  6 


0.4 

0.2 

0 


0.4 

0.2 

0 


0.2  04  0©  o.a 
Mode  7 


0.4 

0.2 

0 

( 

J  0  2  0  4  0.6  0.3 

Mode  S 

0.4 

- 

0.2 

- r * 

n 

0  02  04  0.6  08 


0  0.2  0.4  0.6  0.8 


Figure  7.  8  POD  Modes. 

simulation  at  each  time  step  onto  the  POD  basis.  The 
resulting  first  5  temporal  coefficients  of  the  full  order 
model  are  compared  to  those  predicted  from  the  13-th 
order  reduced  order  model  output. 

In  Figure  8,  dashed  lines  denote  the  linear  combination  of 
POD  modes  restricted  to  the  boundary.  Solid  lines  denote 


the  boundary  inputs  defined  by 

usit)  =  sin  | 

3”) 

(82) 

uj(t)  =  sin 

[lm) 

(83) 

Figure  8.  Boundary  Condition  Accuracy. 


The  results  obtained  for  the  first  five  temporal  coefficients 
are  shown  in  Figure  9.  In  that  Figure,  solid  curves  denote 
values  of  temporal  coefficients  obtained  from  the  projec¬ 
tion.  Dashed  curves  denote  the  output  of  the  reduced  or¬ 
der  model.  As  seen  in  Figure  9,  very  good  agreement  is 
seen  between  the  full  and  reduced  solutions  even  though 
the  open-loop  input  considered  was  not  specifically  incor¬ 
porated  in  the  snapshot  ensemble. 


Figure  9.  Projected  and  POD  Model  Coefficients. 

In  Figure  10,  we  compare  the  full  order  solution  w(t,x)  of 
the  Burgers’  equation  with  the  solution  based  on  the  13-th 
order  model  wr(t,x),  i.e., 


As  can  be  seen  in  Figure  8,  there  is  very  good  agreement  40 

between  the  boundary  conditions  specified  for  the  full  order  Hv(f,x)  =  £y;(r)<Mx)  (84) 

system  and  the  linear  combination  of  POD  modes  restricted  i=  l 

to  the  boundary. 
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equation 


0  0.2  0.4  0.6  0.8 


Figure  10.  Full  and  Reduced  Order  Models’  Re¬ 
sponses. 


Note  that  despite  the  fact  that  only  a  13-th  order  model  is 
used  the  agreement  between  the  two  responses  is  good. 

4.3.  H°°  Control 

In  this  section  we  consider  the  design  of  an  H°°  con¬ 
troller  for  the  tracking  problem  for  (77)-(78).  The  moti¬ 
vation  behind  our  choice  is  that  H°°  controllers  are  robust 
against  unmodeled  or  neglected  dynamics,  and  unknown  or 
unmeasurable  disturbances  [21,  22,  23].  As  in  section  3.6, 
a  fixed  reference  signal  wref  (x)  is  specified  for  the  full  or¬ 
der  system.  Projecting  wref(x)  onto  the  POD  basis  yields 
tracking  coefficients  for  the  reduced  order  model,  denoted 
by  aref.  The  tracking  problem  is  depicted  in  Figure  11, 
where  C  is  the  controller  and  P  the  plant  represented  by 
the  dynamical  equation  (25).  The  computation  of  the  H°°  is 
based  on  the  13-th  order  reduced  model  (81).  From  Figure 
1 1,  for  tracking  purposes  the  controlled  output  is  chosen  to 
be  the  error  signal  e  which  is  defined  to  be  the  difference 
between  the  reference  aref  and  the  actual  output  y(t),  i.e., 

e(t)  :=  aref-y(t)  (85) 

The  dynamics  of  the  reduced  model  for  tracking  control 
represented  in  Figure  1 1  are  then  given  by  the  state  space 


z(t)  =  Az(t)  +B\aref  -\-B2u{t)  (86) 

y(t)  —  Cz(^)  T- Dll  O^ref  T ^12^(0 

e(t)  =  —  Cz{t)  +  D21  ocref  +  D22u(t)  (87) 

The  objective  of  the  H°°  controller  C  is  to  stabilize  the 
closed-loop  system  and  minimize  the  effect  of  aref  on  the 
error  e  by  viewing  aref  as  an  unknown  disturbance  in  L 2 
of  unit  ||  •  1 1 2-norm.  From  Figure  11,  in  terms  of  transfer 
function  matrices  of  P  and  C,  the  transfer  matrix  from  aref 
to  e  is  given  by  the  sensitivity  function  Te0Cref  defined  by 

e  =  (J  +  PC)-'an,f  (88) 

=:  Tearef  (89) 

We  compute  the  worst-case  disturbance  transmission  error 


Figure  11.  Block  Diagram  of  the  Closed-Loop  Sys¬ 
tem. 

due  to  aref ,  i.e., 


sup  ||e||2  (90) 

\\arefh<\ 

which  is  given  by  [21,  22,  23] 

sup  ||e||2  =  ess  sup  o{TeareAj(0 ))  (91) 

II<V/I|2<1  0<ffl<co  V  ' 

=:  ll^are/||oo  (92) 

where  esssup  denotes  the  essential  supremum,  and  cr(-) 
the  maximum  singular  value  of  its  argument. 

The  H°°  control  design  reduces  to  the  following  opti¬ 
mization:  Find  C  such  that  the  closed-loop  system  is 
robustly  stable  and 

p.  ■=  mm\\Tearef\\oo  (93) 

The  solution  of  (93)  is  textbook  material.  There  are 
Riccati-based  and  linear  matrix  inequalities  (LMIs)  based 
techniques  to  solve  (93)  [21,  22,  23,  27].  In  this  work,  we 
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use  the  LMI  approach  because  of  its  numerical  robustness 
and  stability.  The  H°°  problem  (93)  is  directly  optimized 
by  solving  the  following  LMI  problem  [27,  28]: 

Minimize  y  over  R  =  RT  and  S  =  ST  such  that 


NU 

0 


N2i 

0 


where  N\2  and  N2\  denote  the  bases  of  the  null  spaces 
of  (§2,  D\. 2)  and  (— C,  D2\),  respectively,  and  I  is  the 
identity  matrix. 
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Solving  the  LMI  optimization  result  in  an  optimal 
/i  =  1  and  an  optimal  controller  C  of  order  13.  Closing  the 
loop  on  the  full  order  Galerkin  model  using  this  controller 
results  in  the  responses  for  the  temporal  coefficients  a 
shown  in  Figure  12,  where  we  have  plotted  the  first  8  a’s. 
In  that  Figure,  solid  curves  represent  the  values  of  the  8 
temporal  coefficients  a,  and  dashed  lines  represent  the 
reference  aref.  Note  that  excellent  tracking  is  achieved 
with  virtually  zero  steady  state  error. 


Figure  12.  Full  Order  Closed-Loop  System  Tracking 
Response 

quires  linearization  and  is  only  valid  in  the  vicinity  of  the 
equilibrium  for  which  it  was  developed.  Nonetheless,  given 
the  success  using  traditional  POD  and  balanced  truncation 
for  the  Burgers  equation,  both  methods  show  promise  for 
use  with  the  Navier-Stokes  equations. 


5.  Conclusions 

For  the  two-dimensional  Burgers  equation,  traditional 
proper  orthogonal  decomposition  and  balanced  truncation 
using  POD  modes  are  solid  approaches  for  order  reduc¬ 
tion  with  the  goal  of  control  design.  Each  method  has  par¬ 
ticular  strengths  that  make  it  well- suited  for  this  class  of 
problems.  Traditional  POD  maintains  nonlinearities  and 
therefore  allows  the  use  of  the  reduced-order  model  to  as¬ 
sess  the  impact  of  the  nonlinearities  on  the  performance  of 
the  controller  in  reduced-order  simulations.  It  also  permits 
the  use  of  nonlinear  control,  although  that  was  not  a  part 
of  the  current  study.  Balanced  truncation  provides  much 
needed  insight  into  which  modes  should  be  maintained  in 
the  reduced-order  model.  Each  method  also  presents  chal¬ 
lenges.  Traditional  POD,  while  providing  modes  that  are 
physically  meaningful,  does  not  always  provide  modes  that 
are  observable  or  controllable.  Balanced  truncation  re- 


6.  Future  Work 

Application  of  the  order-reduction  methods  will  be  ap¬ 
plied  to  a  two-dimensional,  simply-hinged  airfoil  with  mul¬ 
tiple  synthetic  jets  in  the  vicinity  of  and  on  the  flap  to  elim¬ 
inate  separation  and  generate  high  lift.  A  well-designed 
control  law  is  expected  to  provide  the  correct  phase  of  the 
actuators  with  respect  to  each  other  and  the  flow  field  dy¬ 
namics. 
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